Section 4: snRNAseq analysis

Arun Seetharam, Ha Vu

Tuteja Lab

2021-08-14

snRNAseq analyses

Prerequisites

Packages required for this analyses are loaded as follows:

setwd("~/github/amnion.vs.other_RNASeq")
library(Seurat)
library(SeuratWrappers)
library(knitr)
library(kableExtra)
library(ggplot2)
library(cowplot)
library(patchwork)
library(metap)
library(multtest)
library(gridExtra)
library(dplyr)
library(stringr)
library(TissueEnrich)
library(gprofiler2)
library(tidyverse)
library(enhancedDimPlot)
library(calibrate)
library(ggrepel)
library(dittoSeq)
library(ComplexHeatmap)
library(scales)
library(plotly)

Import datasets

The counts data and its associated metadata (coldata) are imported for analyses

experiment_name = "BAP"
dataset_loc <- "~/TutejaLab/expression-data"
ids <- c("5pcO2_r1", "5pcO2_r2", "20pcO2_r1", "20pcO2_r2", "nCT_D5", "nCT_D10", "nTE_D2", "nTE_D3")
d10x.data <- sapply(ids, function(i){
  d10x <- Read10X(file.path(dataset_loc,i,"filtered_feature_bc_matrix"))
  colnames(d10x) <- paste(sapply(strsplit(colnames(d10x),split="-"), '[[' , 1L ), i, sep="-")
  d10x
})
experiment.data <- do.call("cbind", d10x.data)
bapd8.combined <- CreateSeuratObject(
  experiment.data,
  project = "BAPd8",
  min.cells = 10,
  min.genes = 200,
  names.field = 2,
  names.delim = "\\-")
bapd8.temp <- bapd8.combined
bapd8.combined$log10GenesPerUMI <- log10(bapd8.combined$nFeature_RNA) / log10(bapd8.combined$nCount_RNA)
bapd8.combined$mitoRatio <- PercentageFeatureSet(object = bapd8.combined, pattern = "^MT-")
bapd8.combined$mitoRatio <- bapd8.combined@meta.data$mitoRatio / 100
metadata <- bapd8.combined@meta.data
metadata$cells <- rownames(metadata)
metadata <- metadata %>%
  dplyr::rename(seq_folder = orig.ident,
                nUMI = nCount_RNA,
                nGene = nFeature_RNA, 
                seq_folder = orig.ident)
mt <- ggplot(dat = metadata, aes(x=nUMI, y=nGene, color=mitoRatio)) + 
  geom_point(alpha = 0.5) + 
  scale_colour_gradient(low = "gray90", high = "black") + labs(colour="MT ratio") +
  theme_bw(base_size = 12) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    strip.background = element_blank(),
    panel.border = element_rect(colour = "black")) + 
  xlab("RNA counts") + ylab("Gene counts") +
  stat_smooth(method=lm) +
  facet_wrap(~seq_folder, labeller = labeller(seq_folder = 
                                                c("20pcO2_r1" = "20% Oxygen (rep1)",
                                                  "20pcO2_r2" = "20% Oxygen (rep2)",
                                                  "5pcO2_r1" = "5% Oxygen (rep1)",
                                                  "5pcO2_r2" = "5% Oxygen (rep2)",
                                                  "nCT_D5" = "nCT day 5",
                                                  "nCT_D10" = "nCT day 10",
                                                  "nTE_D2" = "nTE day 2",
                                                  "nTE_D3" = "nTE day 3"))) +
  scale_y_continuous(label=comma) +
  scale_x_continuous(label=comma) 
mt
#> `geom_smooth()` using formula 'y ~ x'
FigX: fig legend

FigX: fig legend

ncells <- ggplot(metadata, aes(x=seq_folder, fill=seq_folder)) + 
  geom_bar() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme(plot.title = element_text(hjust=0.5, face="bold")) +
  ggtitle("Number of Cells")
ncells
FigX: fig legend

FigX: fig legend

density of cells per sample

dcells <- ggplot(metadata, aes(color=seq_folder, x=nUMI, fill= seq_folder)) + 
  geom_density(alpha = 0.2) + 
  scale_x_log10() + 
  theme_classic() +
  ylab("Cell density") +
  geom_vline(xintercept = 500)
dcells
FigX: fig legend

FigX: fig legend

ngenes <- ggplot(metadata, aes(x=seq_folder, y=log10(nGene), fill=seq_folder)) + 
  geom_boxplot() + 
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme(plot.title = element_text(hjust=0.5, face="bold")) +
  ggtitle("NCells vs NGenes")
ngenes
FigX: fig legend

FigX: fig legend

# transcritps
dtranscripts <- ggplot(metadata, aes(x=log10GenesPerUMI, color = seq_folder, fill=seq_folder)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  geom_vline(xintercept = 0.8)
dtranscripts
FigX: fig legend

FigX: fig legend

# mito ratio
mtratio <- ggplot(metadata, aes(color=seq_folder, x=mitoRatio, fill=seq_folder)) + 
  geom_density(alpha = 0.2) + 
  scale_x_log10() + 
  theme_classic() +
  geom_vline(xintercept = 0.2)
mtratio
#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning: Removed 21 rows containing non-finite values (stat_density).
FigX: fig legend

FigX: fig legend

bapd8.combined <- bapd8.temp
df <- bapd8.combined@meta.data
df$replicate <- NA
df$replicate[which(str_detect(df$orig.ident, "5pcO2"))] <- "5pcO2"
df$replicate[which(str_detect(df$orig.ident, "20pcO2"))] <- "20pcO2"
df$replicate[which(str_detect(df$orig.ident, "nCT_"))] <- "nCT"
df$replicate[which(str_detect(df$orig.ident, "nTE_"))] <- "nTE"
bapd8.combined@meta.data <- df
bapd8.combined[["percent.mt"]] <- PercentageFeatureSet(bapd8.combined, pattern = "^MT-")
p <- VlnPlot(bapd8.combined, features = "nFeature_RNA", pt.size = 1) + 
  geom_hline(yintercept=200, color = "red", size=1) +
  geom_hline(yintercept=10000, color = "red", size=1) +
  theme(legend.position = "none")
q <- VlnPlot(bapd8.combined, features = "nCount_RNA", pt.size = 1) +
  theme(legend.position = "none")
r <- VlnPlot(bapd8.combined, features = "percent.mt", pt.size = 1) +
  geom_hline(yintercept=15, color = "red", size=1) +
  theme(legend.position = "none")
panel_plot <- plot_grid(p, q, r, labels=c("A", "B", "C"), ncol=3, nrow = 1)
panel_plot
FigX: fig legend

FigX: fig legend

B1 <- FeatureScatter(bapd8.combined, feature1 = "nCount_RNA", feature2 = "percent.mt")
B2 <- FeatureScatter(bapd8.combined, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
bapd8.combined <- subset(bapd8.combined, subset = nFeature_RNA > 200 & nFeature_RNA < 10000 & percent.mt < 25)
I1 <- FeatureScatter(bapd8.combined, feature1 = "nCount_RNA", feature2 = "percent.mt")
I2 <- FeatureScatter(bapd8.combined, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
bapd8.combined <- subset(bapd8.combined, subset = nFeature_RNA > 200 & nFeature_RNA < 10000 & percent.mt < 15)
A1 <- FeatureScatter(bapd8.combined, feature1 = "nCount_RNA", feature2 = "percent.mt")
A2 <- FeatureScatter(bapd8.combined, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
B <- B1 | B2
I <- I1 | I2
A <- A1 | A2
panel_plot <- plot_grid(B1, B2, I1, I2, A1, A2, labels=c("A", "B", "C", "D", "E", "F"), ncol=2, nrow = 3)
panel_plot
FigX: fig legend

FigX: fig legend

counts <- GetAssayData(object = bapd8.combined, slot = "counts")
counts <- counts[grep(pattern = "^MT", x = rownames(counts), invert = TRUE),]
counts <- counts[grep(pattern = "^MT", x = rownames(counts), invert = TRUE),]
counts <- counts[grep(pattern = "^RPL", x = rownames(counts), invert = TRUE),]
counts <- counts[grep(pattern = "^RPS", x = rownames(counts), invert = TRUE),]
counts <- counts[grep(pattern = "^MRPS", x = rownames(counts), invert = TRUE),]
counts <- counts[grep(pattern = "^MRPL", x = rownames(counts), invert = TRUE),]
keep_genes <- Matrix::rowSums(counts) >= 10
filtered_counts <- counts[keep_genes, ]
bapd8.fcombined <- CreateSeuratObject(filtered_counts, meta.data = bapd8.combined@meta.data)
bapd8.fcombined@meta.data <- bapd8.fcombined@meta.data[1:4]
bapd8.combined <- bapd8.fcombined
bapd8.list <- SplitObject(bapd8.combined, split.by = "orig.ident")
bapd8.list <- lapply(X = bapd8.list, FUN = function(x) {
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
bapd8.anchors <- FindIntegrationAnchors(object.list = bapd8.list, dims = 1:20)
bapd8.integrated <- IntegrateData(anchorset = bapd8.anchors, dims = 1:20)
DefaultAssay(bapd8.integrated) <- "integrated"
bapd8.integrated <- ScaleData(bapd8.integrated, verbose = FALSE)
bapd8.integrated <- RunPCA(bapd8.integrated, npcs = 30, verbose = FALSE)
bapd8.integrated <- RunUMAP(bapd8.integrated, reduction = "pca", dims = 1:20)
#> Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
#> To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
#> This message will be shown once per session
bapd8.integrated <- FindNeighbors(bapd8.integrated, reduction = "pca", dims = 1:20)
bapd8.integrated <- FindClusters(bapd8.integrated, resolution = 0.5)
#> Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
#> 
#> Number of nodes: 17208
#> Number of edges: 610547
#> 
#> Running Louvain algorithm...
#> Maximum modularity in 10 random starts: 0.8716
#> Number of communities: 13
#> Elapsed time: 1 seconds
num.clusters <- nlevels(bapd8.integrated$seurat_clusters)
num.clusters
#> [1] 13
VizDimLoadings(bapd8.integrated, dims = 1:4, reduction = "pca")
FigX: fig legend

FigX: fig legend

bapd8.integrated <- JackStraw(bapd8.integrated, num.replicate = 100)
bapd8.integrated <- ScoreJackStraw(bapd8.integrated, dims = 1:20)
elbow <- ElbowPlot(bapd8.integrated)
jack <- JackStrawPlot(bapd8.integrated, dims = 1:20)
panel_plot <- plot_grid(elbow, jack, labels = c('A', 'B'), ncol = 2)
#> Warning: Removed 28000 rows containing missing values (geom_point).
panel_plot
FigX: fig legend

FigX: fig legend

df <- bapd8.integrated@meta.data
df$new_clusters <- as.factor(as.numeric(df$seurat_clusters))
bapd8.integrated@meta.data <- df
Idents(bapd8.integrated) <- "new_clusters"
A = enhancedDimPlot(object = bapd8.integrated,
                    grouping_var = 'ident',
                    reduction = "umap",
                    label = TRUE,
                    pt.size = 1,
                    alpha = 0.5) +
  ggtitle("A") +
  xlab("UMAP_1") +
  ylab("UMAP_2") +
  theme_classic() + 
  theme(legend.position = "none", 
        plot.title = element_text(face = "bold"))

B <- enhancedDimPlot(object = bapd8.integrated,
                     grouping_var = 'replicate',
                     reduction = "umap",
                     label = FALSE,
                     pt.size = 1,
                     alpha = 0.4) +
  ggtitle("B") +
  xlab("UMAP_1") +
  ylab("UMAP_2") +
  theme_classic() + 
  theme(legend.justification = c(1, 1),
        legend.position = c(1, 1),
        plot.title = element_text(face = "bold")) +
  scale_colour_manual(name = "Conditions", 
                      labels = c(expression(paste('20% ', 'O'[2])),
                                 expression(paste('5% ', 'O'[2])),
                                 'nCT',
                                 'nTE'), 
                      values = c("20pcO2" = "#DA3C96",
                                 "5pcO2" = "#A90065",
                                 "nCT" = "#FFD74D",
                                 "nTE" = "#9BC13C" )) +
  scale_fill_manual(name = "Conditions", 
                    labels = c(expression(paste('20% ', 'O'[2])),
                               expression(paste('5% ', 'O'[2])),
                               'nCT',
                               'nTE'), 
                    values = c("20pcO2" = "#DA3C96",
                                 "5pcO2" = "#A90065",
                                 "nCT" = "#FFD74D",
                                 "nTE" = "#9BC13C")) +
  scale_linetype_manual(values = "blank")

C <- enhancedDimPlot(object = bapd8.integrated,
                     grouping_var = 'orig.ident',
                     reduction = "umap",
                     label = FALSE,
                     pt.size = 1,
                     alpha = 0.4) +
  ggtitle("C") +
  xlab("UMAP_1") +
  ylab("UMAP_2") +
  theme_classic() + 
  theme(legend.justification = c(1, 1),
        legend.position = c(1, 1),
        plot.title = element_text(face = "bold")) +
  scale_colour_manual(name = "Replicates", 
                      labels = c(expression(paste('20% ', 'O'[2], ' rep1')),
                                 expression(paste('20% ', 'O'[2], ' rep2')),
                                 expression(paste('5% ', 'O'[2], ' rep1')),
                                 expression(paste('5% ', 'O'[2], ' rep1')),
                                 "nCT day 5",
                                 "nCT day 10",
                                 "nTE day 3",
                                 "nTE day 5"), 
                      values = c("20pcO2_r1" = "#0571b0",
                                 "20pcO2_r2" = "#92c5de",
                                 "5pcO2_r1" = "#ca0020",
                                 "5pcO2_r2" = "#f4a582",
                                 "nCT_D5" = "#d133ff",
                                 "nCT_D10" = "#ff33f6",
                                 "nTE_D2" = "#33ffa2",
                                 "nTE_D3" = "#5bff33")) +
  scale_fill_manual(name = "Replicates", 
                    labels = c(expression(paste('20% ', 'O'[2], ' rep1')),
                               expression(paste('20% ', 'O'[2], ' rep2')),
                               expression(paste('5% ', 'O'[2], ' rep1')),
                               expression(paste('5% ', 'O'[2], ' rep1')),
                               "nCT day 5",
                               "nCT day 10",
                               "nTE day 3",
                               "nTE day 5"), 
                    values = c("20pcO2_r1" = "#0571b0",
                               "20pcO2_r2" = "#92c5de",
                               "5pcO2_r1" = "#ca0020",
                               "5pcO2_r2" = "#f4a582",
                               "nCT_D5" = "#d133ff",
                               "nCT_D10" = "#ff33f6",
                               "nTE_D2" = "#33ffa2",
                               "nTE_D3" = "#5bff33")) +
  scale_linetype_manual(values = "blank")
panel_plot <- plot_grid(A, B, ncol=2, nrow = 1)
panel_plot
FigX: fig legend

FigX: fig legend

DimPlot(object = bapd8.integrated,
                    split.by = 'orig.ident', ncol = 4) +
  xlab("UMAP_1") +
  ylab("UMAP_2") +
  theme_classic() + 
  theme(legend.position = "none")
FigX: fig legend

FigX: fig legend

A = enhancedDimPlot(object = bapd8.integrated,
                    grouping_var = 'ident',
                    reduction = "umap",
                    label = FALSE,
                    pt.size = 1,
                    alpha = 0.5) +
  xlab("UMAP_1") +
  ylab("UMAP_2") +
  theme_classic() + 
  theme(legend.position = "none")
ggplotly(A)

FigX: fig legend

DefaultAssay(bapd8.integrated) <- "RNA"
for (i in 1:num.clusters){
  try({
    cluster.markers.all <- FindMarkers(bapd8.integrated, ident.1 = i)
    cluster.markers.filtered <- cluster.markers.all %>% filter(avg_log2FC >= 0.584962501) %>% filter(p_val_adj <= 0.05) %>% arrange(desc(avg_log2FC))
    markers.filtered.names <- rownames(cluster.markers.filtered)
    assign( paste("cluster.marker.names", i, sep = "."), markers.filtered.names ) 
  })
}
fullCounts <- tibble(
    cluster = bapd8.integrated@meta.data$new_clusters,
    cell_type = bapd8.integrated@meta.data$orig.ident) %>%
    dplyr::group_by(cluster, cell_type) %>%
    dplyr::count() %>%
    dplyr::group_by(cluster) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(cluster=paste0("Cluster_",cluster))
fullCounts <- fullCounts %>% 
  group_by(cell_type) %>% 
  mutate(cell_type_sum = sum(n)) %>% 
  mutate(percent=(n*100)/cell_type_sum)
list2env(split(fullCounts, fullCounts$cluster), envir = .GlobalEnv)
#> <environment: R_GlobalEnv>
mybarplot <- function(pdata, i) {
  ggplot(data=pdata, aes(x=cell_type, y=percent, fill = cell_type, alpha = 0.5)) +
    geom_bar(stat="identity") +
    theme_classic(base_size = 12) +
    theme(legend.position="none",
          axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    ggtitle(paste("Cluster", i, "cell composition")) + ylab("% cells in cluster") + xlab("") +
    scale_colour_manual(name = "Replicates", 
                      labels = c(expression(paste('20% ', 'O'[2], ' rep1')),
                                 expression(paste('20% ', 'O'[2], ' rep2')),
                                 expression(paste('5% ', 'O'[2], ' rep1')),
                                 expression(paste('5% ', 'O'[2], ' rep1')),
                                 "nCT day 5",
                                 "nCT day 10",
                                 "nTE day 3",
                                 "nTE day 5"), 
                      values = c("20pcO2_r1" = "#DA3C96",
                               "20pcO2_r2" = "#DA3C96",
                               "5pcO2_r1" = "#A90065",
                               "5pcO2_r2" = "#A90065",
                               "nCT_D5" = "#FFD74D",
                               "nCT_D10" = "#FFD74D",
                               "nTE_D2" = "#9BC13C",
                               "nTE_D3" = "#9BC13C")) +
  scale_fill_manual(name = "Replicates", 
                    labels = c(expression(paste('20% ', 'O'[2], ' rep1')),
                               expression(paste('20% ', 'O'[2], ' rep2')),
                               expression(paste('5% ', 'O'[2], ' rep1')),
                               expression(paste('5% ', 'O'[2], ' rep1')),
                               "nCT day 5",
                               "nCT day 10",
                               "nTE day 3",
                               "nTE day 5"), 
                    values = c("20pcO2_r1" = "#DA3C96",
                               "20pcO2_r2" = "#DA3C96",
                               "5pcO2_r1" = "#A90065",
                               "5pcO2_r2" = "#A90065",
                               "nCT_D5" = "#FFD74D",
                               "nCT_D10" = "#FFD74D",
                               "nTE_D2" = "#9BC13C",
                               "nTE_D3" = "#9BC13C")) +
  scale_linetype_manual(values = "blank")
}
# colors for each clusters defined
ggplotColours <- function(n = 6, h = c(0, 360) + 15){
  if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
  hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}
color_list <- ggplotColours(n=13)
grouped_violinPlots <- function(markersfile, clusternumber, seuratobject = bapd8.integrated) {
  dittoPlotVarsAcrossGroups(seuratobject, markersfile, 
                            group.by = "new_clusters", main = paste("Cluster ", clusternumber, " markers expression"), 
                            xlab = "", 
                            ylab = "Mean z-score expression", 
                            x.labels = c("Cluster 1", "Cluster 2", "Cluster 3", "Cluster 4",
                                         "Cluster 5", "Cluster 6", "Cluster 7", "Cluster 8",
                                         "Cluster 9", "Cluster 10", "Cluster 11", "Cluster 12",
                                         "Cluster 13"), 
                            vlnplot.lineweight = 0.5, 
                            legend.show = FALSE, 
                            jitter.size = 0.5, 
                            color.panel = color_list)
}
l <- load(file = "~/TutejaLab/PlacentaEnrich/combine-test-expression1.Rdata")
humanGeneMapping <- dataset$GRCH38$humanGeneMapping
d <- dataset$PlacentaDeciduaBloodData
data <- d$expressionData
cellDetails <- d$cellDetails
runpce <- function(inputgenelist, clusternumber, shade) {
  inputGenes<-toupper(inputgenelist)
  humanGene<-humanGeneMapping[humanGeneMapping$Gene.name %in% inputGenes,]
  inputGenes<-humanGene$Gene
  expressionData<-data[intersect(row.names(data),humanGeneMapping$Gene),]
  se<-SummarizedExperiment(assays = SimpleList(as.matrix(expressionData)),rowData = row.names(expressionData),colData = colnames(expressionData))
  cellSpecificGenesExp<-teGeneRetrieval(se,expressedGeneThreshold = 1)
  print(length(inputGenes))
  gs<-GeneSet(geneIds=toupper(inputGenes))
  output2<-teEnrichmentCustom(gs,cellSpecificGenesExp)
  enrichmentOutput<-setNames(data.frame(assay(output2[[1]]),row.names = rowData(output2[[1]])[,1]),colData(output2[[1]])[,1])
  row.names(cellDetails)<-cellDetails$RName
  enrichmentOutput$Tissue<- cellDetails[row.names(enrichmentOutput),"CellName"]
  ggplot(data = enrichmentOutput, mapping = aes(x = reorder (Tissue, -Log10PValue), Log10PValue)) +
    geom_bar(stat = "identity", color = shade, fill = shade) + theme_classic(base_size = 10) + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size = 10), plot.margin = unit(c(1,1,1,3), "cm")) +
    ggtitle(paste0("Cluster ", clusternumber, " PCE")) +
    ylab("-log10 p-value (adj.)") +
    xlab("") +
    scale_y_continuous(expand = expansion(mult = c(0, .1)))
}
pce <- runpce(cluster.marker.names.1, 1, color_list[1])
#> [1] 68
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_1, 1)
violin <- grouped_violinPlots(cluster.marker.names.1, 1)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.2, 2, color_list[2])
#> [1] 44
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_2, 2)
violin <- grouped_violinPlots(cluster.marker.names.2, 2)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.3, 3, color_list[3])
#> [1] 244
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_3, 3)
violin <- grouped_violinPlots(cluster.marker.names.3, 3)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.4, 4, color_list[4])
#> [1] 210
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_4, 4)
violin <- grouped_violinPlots(cluster.marker.names.4, 4)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.5, 5, color_list[5])
#> [1] 32
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_5, 5)
violin <- grouped_violinPlots(cluster.marker.names.5, 5)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.6, 6, color_list[6])
#> [1] 355
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_6, 6)
violin <- grouped_violinPlots(cluster.marker.names.6, 6)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.7, 7, color_list[7])
#> [1] 38
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_7, 7)
violin <- grouped_violinPlots(cluster.marker.names.7, 7)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.1, 8, color_list[8])
#> [1] 68
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_8, 8)
violin <- grouped_violinPlots(cluster.marker.names.8, 8)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.1, 9, color_list[9])
#> [1] 68
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_9, 9)
violin <- grouped_violinPlots(cluster.marker.names.9, 9)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.1, 10, color_list[10])
#> [1] 68
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_10, 10)
violin <- grouped_violinPlots(cluster.marker.names.10, 10)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.1, 11, color_list[11])
#> [1] 68
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_11, 11)
violin <- grouped_violinPlots(cluster.marker.names.11, 11)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.1, 12, color_list[12])
#> [1] 68
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_12, 12)
violin <- grouped_violinPlots(cluster.marker.names.12, 12)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

pce <- runpce(cluster.marker.names.1, 13, color_list[13])
#> [1] 68
#> No background list provided. Using all the
#>                 genes as background.
count <- mybarplot(Cluster_13, 13)
violin <- grouped_violinPlots(cluster.marker.names.13, 13)
toprow <- plot_grid(count, violin, labels = c('A', 'B'), align = 'h')
panel_plot <- plot_grid(toprow, pce, labels = c('', 'C'), ncol = 1, rel_heights = c(1, 1.3))
panel_plot
FigX: fig legend

FigX: fig legend

sessionInfo()
#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19042)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.1252 
#> [2] LC_CTYPE=English_United States.1252   
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.1252    
#> 
#> attached base packages:
#>  [1] grid      stats4    parallel  stats     graphics  grDevices utils    
#>  [8] datasets  methods   base     
#> 
#> other attached packages:
#>  [1] plotly_4.9.3                scales_1.1.1               
#>  [3] ComplexHeatmap_2.6.2        dittoSeq_1.2.6             
#>  [5] ggrepel_0.9.1               calibrate_1.7.7            
#>  [7] MASS_7.3-54                 enhancedDimPlot_0.0.0.9100 
#>  [9] forcats_0.5.1               purrr_0.3.4                
#> [11] readr_1.4.0                 tidyr_1.1.3                
#> [13] tibble_3.1.1                tidyverse_1.3.1            
#> [15] gprofiler2_0.2.0            TissueEnrich_1.10.1        
#> [17] GSEABase_1.52.1             graph_1.68.0               
#> [19] annotate_1.68.0             XML_3.99-0.6               
#> [21] AnnotationDbi_1.52.0        SummarizedExperiment_1.20.0
#> [23] GenomicRanges_1.42.0        GenomeInfoDb_1.26.7        
#> [25] IRanges_2.24.1              S4Vectors_0.28.1           
#> [27] MatrixGenerics_1.2.1        matrixStats_0.58.0         
#> [29] ensurer_1.1                 stringr_1.4.0              
#> [31] dplyr_1.0.7                 gridExtra_2.3              
#> [33] multtest_2.46.0             Biobase_2.50.0             
#> [35] BiocGenerics_0.36.1         metap_1.4                  
#> [37] patchwork_1.1.1             cowplot_1.1.1              
#> [39] ggplot2_3.3.5               kableExtra_1.3.4           
#> [41] knitr_1.33                  SeuratWrappers_0.3.0       
#> [43] SeuratObject_4.0.1          Seurat_4.0.2               
#> 
#> loaded via a namespace (and not attached):
#>   [1] scattermore_0.7             bit64_4.0.5                
#>   [3] irlba_2.3.3                 multcomp_1.4-17            
#>   [5] DelayedArray_0.16.3         data.table_1.14.0          
#>   [7] rpart_4.1-15                RCurl_1.98-1.3             
#>   [9] generics_0.1.0              TH.data_1.0-10             
#>  [11] RSQLite_2.2.7               RANN_2.6.1                 
#>  [13] future_1.21.0               bit_4.0.4                  
#>  [15] mutoss_0.1-12               spatstat.data_2.1-0        
#>  [17] webshot_0.5.2               xml2_1.3.2                 
#>  [19] lubridate_1.7.10            httpuv_1.6.1               
#>  [21] assertthat_0.2.1            xfun_0.22                  
#>  [23] hms_1.1.0                   jquerylib_0.1.4            
#>  [25] evaluate_0.14               promises_1.2.0.1           
#>  [27] fansi_0.4.2                 dbplyr_2.1.1               
#>  [29] readxl_1.3.1                igraph_1.2.6               
#>  [31] DBI_1.1.1                   tmvnsim_1.0-2              
#>  [33] htmlwidgets_1.5.3           spatstat.geom_2.1-0        
#>  [35] paletteer_1.3.0             ellipsis_0.3.2             
#>  [37] crosstalk_1.1.1             RSpectra_0.16-0            
#>  [39] backports_1.2.1             bookdown_0.22              
#>  [41] deldir_0.2-10               vctrs_0.3.8                
#>  [43] SingleCellExperiment_1.12.0 remotes_2.3.0              
#>  [45] Cairo_1.5-12.2              ROCR_1.0-11                
#>  [47] abind_1.4-5                 cachem_1.0.5               
#>  [49] withr_2.4.2                 sctransform_0.3.2          
#>  [51] goftest_1.2-2               mnormt_2.0.2               
#>  [53] svglite_2.0.0               cluster_2.1.2              
#>  [55] lazyeval_0.2.2              crayon_1.4.1               
#>  [57] labeling_0.4.2              edgeR_3.32.1               
#>  [59] pkgconfig_2.0.3             nlme_3.1-152               
#>  [61] rlang_0.4.11                globals_0.14.0             
#>  [63] lifecycle_1.0.0             miniUI_0.1.1.1             
#>  [65] sandwich_3.0-1              mathjaxr_1.4-0             
#>  [67] modelr_0.1.8                rsvd_1.0.5                 
#>  [69] cellranger_1.1.0            polyclip_1.10-0            
#>  [71] lmtest_0.9-38               Matrix_1.3-3               
#>  [73] zoo_1.8-9                   reprex_2.0.0               
#>  [75] ggridges_0.5.3              GlobalOptions_0.1.2        
#>  [77] pheatmap_1.0.12             png_0.1-7                  
#>  [79] viridisLite_0.4.0           rjson_0.2.20               
#>  [81] bitops_1.0-7                KernSmooth_2.23-20         
#>  [83] blob_1.2.1                  shape_1.4.6                
#>  [85] parallelly_1.25.0           memoise_2.0.0              
#>  [87] magrittr_2.0.1              plyr_1.8.6                 
#>  [89] ica_1.0-2                   zlibbioc_1.36.0            
#>  [91] compiler_4.0.5              RColorBrewer_1.1-2         
#>  [93] clue_0.3-59                 plotrix_3.8-1              
#>  [95] fitdistrplus_1.1-5          cli_2.5.0                  
#>  [97] XVector_0.30.0              listenv_0.8.0              
#>  [99] pbapply_1.4-3               mgcv_1.8-35                
#> [101] tidyselect_1.1.1            stringi_1.6.2              
#> [103] highr_0.9                   yaml_2.2.1                 
#> [105] locfit_1.5-9.4              sass_0.4.0                 
#> [107] tools_4.0.5                 future.apply_1.7.0         
#> [109] circlize_0.4.12             rstudioapi_0.13            
#> [111] farver_2.1.0                rmdformats_1.0.3           
#> [113] Rtsne_0.15                  digest_0.6.27              
#> [115] BiocManager_1.30.15         shiny_1.6.0                
#> [117] Rcpp_1.0.6                  broom_0.7.6                
#> [119] later_1.2.0                 RcppAnnoy_0.0.18           
#> [121] httr_1.4.2                  Rdpack_2.1.1               
#> [123] colorspace_2.0-1            rvest_1.0.0                
#> [125] fs_1.5.0                    tensor_1.5                 
#> [127] reticulate_1.20             splines_4.0.5              
#> [129] uwot_0.1.10                 rematch2_2.1.2             
#> [131] sn_2.0.0                    spatstat.utils_2.1-0       
#> [133] systemfonts_1.0.2           xtable_1.8-4               
#> [135] jsonlite_1.7.2              R6_2.5.0                   
#> [137] TFisher_0.2.0               pillar_1.6.1               
#> [139] htmltools_0.5.1.1           mime_0.10                  
#> [141] glue_1.4.2                  fastmap_1.1.0              
#> [143] codetools_0.2-18            mvtnorm_1.1-1              
#> [145] utf8_1.2.1                  lattice_0.20-44            
#> [147] bslib_0.2.5.1               spatstat.sparse_2.0-0      
#> [149] numDeriv_2016.8-1.1         leiden_0.3.8               
#> [151] survival_3.2-11             limma_3.46.0               
#> [153] rmarkdown_2.10.1            munsell_0.5.0              
#> [155] GetoptLong_1.0.5            GenomeInfoDbData_1.2.4     
#> [157] haven_2.4.1                 reshape2_1.4.4             
#> [159] gtable_0.3.0                rbibutils_2.1.1            
#> [161] spatstat.core_2.1-2